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For networks of pulse-coupled oscillators with complex connectivity, we demonstrate that in the presence 
of coupling heterogeneity precisely timed periodic firing patterns replace the state of global synchrony that 
exists in homogenous networks only. With increasing disorder, these patterns persist until they reach a critical 
temporal extent that is of the order of the interaction delay. For stronger disorder these patterns cease to exist and 
only asynchronous, aperiodic states are observed. We derive self-consistency equations to predict the precise 
temporal structure of a pattern from the network heterogeneity. Moreover, we show how to design heterogenous 
coupling architectures to create an arbitrary prescribed pattern. 
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Understanding how the structure of a complex network (3] 
determines its dynamics is currently in the focus of research 
in physics, biology and technology |2]. Pulse-coupled os- 
cillators provide a paradigmatic class of models to describe 
a variety of networks that occur in nature, such as popula- 
tions of fireflies, pacemakers cells of the heart, earthquakes, 
or neural networks yl|4[. Synchronization is one of the most 
prevalent kinds of collective dynamics in such networks @1 ■ 
Recent theoretical studies, which analyze conditions for the 
existence and stability of synchronous states, have focused 
on homogenous networks with simple topologies, e.g. global 
couplings or regular lattices |7|]. In nature, however, in- 
tricately structured and heterogenous interactions are ubiqui- 
tous. Previously, aspects of heterogeneity have been studied 
mostly in globally coupled networks liil^l Holl . However, for 
networks with structured connectivity only a few studies exist 
and it is still an open question, how heterogeneity in- 
fluences the dynamics, in particular synchronization, in such 
networks. 

In this Letter we present an exact analysis of the dynam- 
ics of complex networks of pulse-coupled oscillators in the 
presence of coupling heterogeneity. We demonstrate that the 
synchronous state, that exists in homogenous networks, is re- 
placed by precisely timed periodic firing patterns. The tempo- 
ral extent of these patterns grows with the degree of disorder. 
Patterns persist below a critical strength of the disorder be- 
yond which only asynchronous, aperiodic states are observed. 
We show how this transition is controlled by the interaction 
delay. Simple criteria for the stability of firing patterns are 
derived. Furthermore, an approach is presented to predict the 
relative timings of firing events from self-consistency condi- 
tions for the phases in a given network. Conversely, any pre- 
scribed pattern can be created by designing a heterogenous 
coupling architecture in a network of specified connectivity. 

Consider a system of N pulse-coupled oscillators that in- 
teract via directed connections. A matrix C defines the con- 
nectivity of the network, where Cy = 1 if a connection from 
oscillator j to i exists, and CV, = otherwise. The number of 
inputs ki := J2j C%j to everv oscillator i is non-zero, fc; > 1, 
and no further restriction on the network topology is imposed. 
In simulations we consider random graphs in which every di- 
rected connection is present with probability p. 
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FIG. 1: Heterogeneities in a random network (N = 100, p = 0.4, 
b = 2, e = —0.1, r = 0.1) induce desynchronization followed by 
a transition to an asynchronous, aperiodic state, (a) Relative phases 
(ten selected phases shown in black) for different disorder strengths 
Ae (fixed J), (b) Increase of the extent A</> of the pattern shown in 
(a), until the pattern disappears at A<f> c m r (dashed line at A</> = t). 
Inset: Onset of asynchronous state (A(f> ~ 1) for large disorder. 

(c) ,(d) Instantaneous network rate of (c) a periodic (Ae = 0.10) and 

(d) an aperiodic (Ae = 0.15) state in real time computed using a 
triangular sliding time window. Grey histograms on the right of (c) 
and (d) indicate rate distribution (excluding peak at zero network rate 
in (c)). 



The state of an individual oscillator j is represented by a 
phase variable <j>j that increases uniformly in time @], 



d<j)j/dt = 1. 



(1) 



Upon crossing a firing threshold, 4>j{tj) > 1, at time tf an 
oscillator is instantaneously reset to zero, 4»j{t~j) — 0, and a 
pulse is sent. After a delay time r, this pulse is received by all 
oscillators i connected to j (for which Cy = 1) and induces 
an instantaneous phase jump 



H(t f + r)+) = U-'iUiMtf + r)+ 



(2) 
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such that in general the total input is different for every oscil- 
lator. The matrix J specifies the structure of the heterogeneity 
where each element Jy is drawn randomly from the uniform 
distribution on [—1,1]. The degree of heterogeneity is quan- 
tified by the disorder strength Ae. In numerical simulations, 
we take the homogenous part of the coupling strength to be 

The synchronous state © of homogenous networks (Ae = 
0) is replaced by periodic firing patterns in the presence of het- 
erogeneity (Ae > 0), cf. Fig.n After a discontinuous transi- 
tion at a certain critical disorder Ae c these patterns disappear 
and aperiodic, asynchronous states are observed fl5ll . 

The firing patterns emerging for Ae < Ae c are confined to 
a subinterval of the phase axis. The extent of a pattern 
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FIG. 2: The critical extent of a pattern A<j> c as a function of the 
interaction delay r (N — 100, b — 2, £ = —0.1). Each data 
point is obtained for a different heterogeneity J in fully connected 
networks (□) and random networks (O) of different connectivity C 
with p = 0.2. For the random network, data are averaged over 20 
trials (dashed line: fit). Grey area indicates the theoretical second- 
order prediction of the transition zone between A</>f lob ' and A<#™ x . 
Inset: Probability that a periodic pattern emerges from synchronous 
initial conditions for given disorder Ae (histogram of 100 random 
heterogeneities and connectivities, p = 0.2, r = 0.1). 



Here, e^ are the coupling strengths from j to i, which are 
taken to be either purely inhibitory (all e^ < 0) or purely 
excitatory (all e%j > 0). The interaction function U is mono- 
tonically increasing, U' > 0, and represents the subthresh- 
old dynamics of individual oscillators. We consider func- 
tions with a curvature of constant sign, i.e. U"{<j>) > or 
U"(<j>) < for all cf>. This choice characterizes a wide range 
of biologically relevant behaviors, such as neural membrane 
potential dynamics. The qualitative behavior of systems with 
a more complicated interaction function may be derived from 
the results presented here (cf. fl3lo . In simulations, we use 
U (<f>) = ln(l + (e b — l)4>)/b, where b parameterizes the cur- 
vature of U. 

We briefly consider homogenous networks with individual 
coupling strengths eij = iij that are normalized 



N 



(3) 



such that every oscillator i receives the same total input. Such 
networks exhibit a synchronous state, defined by 



(4) 



for all i. The stability of this state has been determined pre- 
viously fl4ll and depends on the interplay between the sign of 
the coupling (excitatory or inhibitory) and the local sign of the 
curvature of the interaction function U fTjll . 

Heterogeneity of the couplings is introduced through 



(taking into account wrap-around effects at the threshold) in- 
creases as the disorder strength Ae is increased, whereas the 
firing order is determined by the structure J of the hetero- 
geneity. If Ae > Ae c the periodic pattern is replaced either 
directly by an asynchronous state (Fig. or it reaches this 
state through a sequence of states consisting of several sepa- 
rated clusters of firing patterns (not shown). 

To understand the origin of this transition, we note that the 
critical disorder Ae c is directly related to the critical pattern 
extent A(j) c before breakdown (cf. Fig.Qb)). A periodic pat- 
tern is always observed for Acj) < r. For globally connected 
networks, the transition to an asynchronous state occurs at 

A(t)f ob - = t (7) 

(see Fig.|2ji. If the system is initialized with a pattern of ex- 
tent A<f) > A0| lob -, some oscillators, which we term critical, 
send pulses that affect oscillators that have not fired within 
the period considered, causing divergence from the periodic 
state. In networks that are not globally connected, the firing 
patterns may persist up to a critical extent Atp c > A<fif loh - 
that depends on the specific network structure. For random 
networks, we estimate an upper bound A(/)™ ax for A(j> c by as- 
suming that the firing times (i.e. the phases) within the pattern 
are uniformlyand independently distributed in an interval of 
length Acj) 1 16]. We first calculate the probability 



P' = l- 



A<p 



1 



k+l 



(8) 



that a given oscillator i is critical, i.e. that of the approximately 
k = pN oscillators j receiving input from i, at least one satis- 
fies <f>i — 4>j > t. In a second step, we calculate the probability 
P = 1 — (1 — P') that at least one of the N oscillators is 
critical. Setting P = 1 guarantees the existence of a criti- 
cal oscillator and results in an implicit formula for A</>™ ax in 
terms of N, p and r. To second order in A</> — r we obtain 



A0™ ax = (1 



(9) 
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as an approximate upper bound on the critical pattern extent 
Acf) c . Numerical data are in good agreement with our theoret- 
ical predictions A0f lob - and A(#° ax (cf. Fig.0. 
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To examine the stability of periodic firing patterns in the 
presence of heterogenous couplings e^, we consider a pattern 



(f>i(zT) = O + &(f>i,o, 



z e 



(10) 



of period T where O represents a common reference phase 
and A0i o defines the relative phase shift of oscillator i within 
the pattern. For simplicity, we order the firing events ac- 
cording to A^i.i > Acf>i t 2 > . . . > A<f>i :ki where A0i.„ 
denotes the phase shift of the oscillator from which oscilla- 
tor i receives its n-th signal within a given period. For ex- 
ample, if f fires first of all oscillators connected to i, then 
A(f> iA = A^v i0 , cf. Eq. GO). 

By definition, all phases increase uniformly in time except 
for two kinds of discrete events: sending and receiving pulses. 
We follow the dynamics of a periodic state, 0, (0) = 0i (T) 
for all i, event by event (cf. UM) starting with a reference 
phase cf>o chosen such that all oscillators have fired but not 
yet received the generated pulses within a given period. The 
collective period of the unperturbed firing pattern is given by 

Acj) 
T i,ki ' 



T = t + 1 + A0j,o - A0 



i.k, 



(ID 



where af := t and 



af n := U l {U{al n _ l + X i,n-i - Xi,n) + e<,n) (12) 

(here with x = A0) are the recursively defined phases of 
oscillator i right after having received the first n of its fc j inputs 
(e» i, ■ ■ ■ , £i.ki, ordered accordingly) in a given period. 
Adding small perturbations i5^o(0) =: 5j,o satisfying 



max \8in\ <-( min |A0 iO - A0,-,o 

i Z \ %,3 



(13) 



to the phases of the oscillators i before firing ensures preserva- 
tion of the ordering of the firing events. Thus, denoting initial 
phases as 0i(O) = 0o + Aipi^, where Ay^o := A0i j o + <5,:.o, 
and sorting these relative phase shifts as before, At^.i > 
Aip it2 > ••• > Aip itki , it follows from A<p iik < Afaj 



that A(p iik < A(p it i implying Atpi 



A0 V 



for 



all i and for all n. Following the evolution of the perturbed 
state event by event we obtain an expression for the time 
T-i = j + 1 — A(fi iki — a^lf to the next firing of i. This 
leads to a nonlinear period-T map 

Si,o(T) =T-(Ti+J + A0 iiO ) =8 iM -a, 



T iZ 



(14) 



of the perturbations, where we used definition d!2t with \ = 
A0 and \ = Ay. Approximating 
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T i,kt 



ki 
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n=l 



(15) 



to first order in 5 ifi , where p 4j „ := Y^iLn {u' {af^ + Afaj - 
A<t>ij +1 ) /U' {aff +1 ) \ forn < fcj andp ijfci = 1, yields 



N 
.7 = 1 



(16) 



Here, the matrix My has diagonal elements Ma = pi$, and 
the non-zero, off-diagonal elements are My = pi >n — Pi t „-i, 
where j is the index of the oscillator sending the n-th signal 
to i (i.e. A<j)jfi = A(j>i, n ). Time translation invariance implies 
£\ = 1 for all i. In addition, if either U" > and e > 
or U" < and i < 0, then p.j j < pj im < 1 for I < m which 
guarantees that all Mij > 0. Using these properties of M it is 
straightforward to show that max^ |<5i,o(^) I < max^ |5i,o(0)|, 
which implies stability of the firing pattern. 

Numerical simulations indicate that firing patterns are de- 
termined by the network topology and the structure of the het- 
erogeneity. How does heterogeneity in a complex network 
control the precise timing of pulses that constitute a pattern? 

In the following we analytically predict the pattern result- 
ing from the underlying coupling architecture. Assume that 
a given coupling architecture induces a pattern dlOt of pe- 
riod T. A change in the coupling matrix e\j = £y +fey below 
the critical disorder leads to a new periodic state ip'^zT') = 
0i(O) + 6(f>ifi of similar period T" that is characterized by the 
phase shift 50 = (<50i,o, ■ ■ ■ , S4>n,o)- Defining A<I>i : „ := 
A0i ; „ + Sipi^n, we pick an arbitrary reference oscillator Z and 
adjust time such that its reference phase is A$i.o = 0. The 
period of oscillator Z in the network defined by the couplings 
e'y is then given by T z = 1+r- Acj) Lkl -5cj) Lkl -er ; A fe * (where 
S<j>ij is an ordering of 50,,o i n analogy to the ordering Afaj 
of A0j.o). After the time T; the phase of an oscillator i ^ Z is 
given by 

ffi) = - + A4>i, ki - A<f> lA + 54>i, ki - S<fn, kl - 

(17) 

Periodicity of this state requires 



Ri(Sc/),Ss) :=#(T,)-#(0) = 



(18) 



for all i, which is an exact, implicit system of algebraic equa- 
tions for 8<p. The collective period T" of the pattern is then 
given by 7} evaluated at the actual phase shifts 5cf>. Expanding 
Ri and T/ to first order in Seij and 6(/>i t o yields 



E 



dR,. 



dS4>jfi 



(0,0) 
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OR, 



dSsi 



(0,0) 



Set. 



(19) 



where i?;(0, 0) = 4 (T) - 0,(0) = by definition. Thus <50 
can be approximated to first order by solving a system of N 
linear equations. This first-order prediction yields a good ap- 
proximation of the actual firing pattern for small disorder Ae 
(see, e.g. Fig. |3ja),(b)) and reasonably describes the qualita- 
tive relative firing times up to Ae = Ae c (cf. Fig.|3jb)). 

The above method may be reversed to design heterogenous 
coupling architectures in order to create a prescribed pattern 
defined by dcf>. Here, we fix the period T' via the reference 
oscillator Z, and solve each of the plane equations i derived by 
linearizing i ll 8t with T' = Ti for a suitable set of Ssij . Two 
examples of networks designed to exhibit the same prescribed 
pattern are shown in Fig.|5Jc). 

In summary, we have demonstrated that for small coupling 
disorder the synchronous state is replaced by precisely timed 
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FIG. 3: Predicting firing patterns and designing coupling architec- 
tures to create a prescribed pattern (6 = 2, e = —0.1, r = 0.1, pre- 
dictions to first order). (a),(b) Comparison of prediction (grey lines) 
and simulation (markers) of selected phases of a pattern (N — 100, 
p — 0.4). (a) Phases shown period by period for Ae = 0.25. Het- 
erogeneity is instantaneously introduced at time t = 0. The attractor 
is reached asymptotically, (b) Pattern for increasing disorder Ae (J 
fixed). Inset: standard deviation of prediction from simulation, (c) 
Coupling strengths (indicated by shading) in two different heteroge- 
nous networks (N = 11, left: p — 0.6, right: p — 0.3) designed 
to create a prescribed firing pattern (middle row, grey). Patterns in 
black are obtained by simulating these designed networks (initialized 
to the synchronous state). 

periodic firing patterns. We have shown how to predict the 
precise timing of pulses from a given coupling architecture, 
and reversely how to design networks in order to create pre- 
scribed patterns. There is a critical disorder strength, which 
we relate to the interaction delay r, beyond which only asyn- 
chronous, aperiodic states are observed. 

Similar behavior is expected for other sources of hetero- 
geneity, such as delay distributions, oscillator specific in- 
teractions Ui(4>) or heterogenous frequencies. The continu- 
ous transition found for small disorder is very different from 
that found previously for threshold-induced synchronization 
where oscillators are split into one synchronous and one asyn- 
chronous subpopulation due to heterogeneity 0|. It has been 
shown recently that periodic firing patterns can be obtained 
and even learned in networks with global inhibition and no 
interaction delays as a perturbation of the asynchronous state 
Jl7ll . In contradistinction, the patterns discussed in this Let- 



ter are not related to (but coexist with) asynchronous states, 
and emerge from the synchronous state due to heterogene- 
ity. The simplicity of the original synchronous state aided us 
to clarify how heterogeneity of networks with a complicated 
structure controls their precise dynamics. These results are 
a first step towards understanding and designing the dynam- 
ics of real world networks that exhibit complex structure and 
heterogeneity. 
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